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ABSTRACT 



We describe analytical and numerical coUisional evolution calculations for 



*^ ' the size distribution of icy bodies in the Kuiper Belt. For a wide range of bulk 

Q ■ properties, initial masses, and orbital parameters, our results yield power-law 

^ ■ cumulative size distributions, A'^^ oc r~'^, with q^ ^ 3.5 for large bodies with 



radii, r > 10-100 km, and Qs ~ 2.5-3 for small bodies with radii, r < 0.1-1 km. 



Qh' The transition between the two power laws occurs at a break radius, r^ ^ 1-30 

O ' km. The break radius is more sensitive to the initial mass in the Kuiper Belt and 

-^ ' the amount of stirring by Neptune than the bulk properties of individual Kuiper 

J^ ■ Belt objects (KBOs). Comparisons with observations indicate that most models 

can explain the observed sky surface density a{m) of KBOs for red magnitudes 

^ ■ i? ~ 22-27. For R < 22 and i? > 28, the model a{m) is sensitive to the amount 

of stirring by Neptune, suggesting that the size distribution of icy planets in the 

outer solar system provides independent constraints on the formation of Neptune. 



Subject headings: planetary systems - solar system: formation - stars: formation 
- circumstellar matter 



1. INTRODUCTION 

The Kuiper Belt is a vast swarm of icy bodies beyond the orbit of Neptune in our solar 
system. Following the discovery of the first Kuiper Belt objects (KBOs) in 1930 (Pluto; 
Tombaugh 1946) and 1992 (1992 QBi; Jewitt & Luu 1993), several groups began large- 
scale surveys to characterize the limits of the Kuiper Belt (e.g., Luu et al 1997; Allen et al. 
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2001; Gladinan et al. 2001; Larsen et al. 2001, and references therein). Today, there are 
800-1000 known KBOs with radii r > 50 km in orbits that extend from 35 AU out to at 
least 150 AU (Luu & Jewitt 2002; Bernstein et al. 2003). The vertical scale height of the 
KBO population is ~ 20-30 deg. The total mass is ~ 0.01-0.1 M® (Luu & Jewitt 2002). 

Observations place numerous constraints on the apparent size distribution of KBOs. 
Deep imaging surveys at i? ~ 22-28 suggest a power-law cumulative size distribution, Nc ex: 
r-'i% with Qo = 3.0 ± 0.5 (Trujillo, Jewitt, & Luu 2001; Luu & Jewitt 2002). Data from 
the ACS on HST suggest a change in the slope of the size distribution at a break radius 
Tb ~ 10-30 km (Bernstein et al. 2003). Dynamical considerations derived from the orbits 
of Pluto-Charon and Jupiter family comets suggest r^ ~ 1-10 km (Duncan et al. 1988; 
Levison & Duncan 1994; Levison & Stern 1995; Duncan et al. 1995; Duncan & Levison 
1997; Ip & Fernandez 1997; Stern, Bottke & Levison 2003). Observations of the optical 
and far-infrared background light require Qo ^ 2.5 for small objects with radii < 0.1-1 
km (Backman, Dasgupta, & Stencel 1995; Stern 1996b; Teplitz et al. 1999; Kenyon & 
Windhorst 2001), suggesting rb > 0.1-1 km. 

These observations provide interesting tests of planet formation theories. In the plan- 
etesimal hypothesis, planetesimals with radii < 1-10 km collide, merge, and grow into larger 
objects. This accretion process yields a power-law size distribution, Nc oc r""^, with qi ~ 
2.8-3.5 for 10-100 km and larger objects (Kenyon & Luu 1999b; Kenyon 2002). As the 
largest objects grow, their gravity stirs up the orbits of leftover planetesimals to the disrup- 
tion velocity. Collisions between leftover planetesimals then produce fragments instead of 
mergers. For objects with radii of 1-10 km and smaller, this process yields a power-law size 
distribution with a shallower slope, qs ~ 2.5 (Stern 1996a; Davis & Farinella 1997; Stern 
& Colwell 1997a,b; Kenyon & Luu 1999a; Kenyon & Bromley 2004). 

Despite uncertainties in the observations and the theoretical calculations, the predicted 
power-law slopes for large and small KBOs agree remarkably well with the data. However, 
many issues remain. The observed sky surface density of the largest objects with radii 
of 300-1000 km and the location of the break in the size distribution are uncertain (see, 
for example, Luu & Jewitt 2002; Bernstein et al. 2003). The data also suggest that 
different dynamical classes of KBOs have different size distributions (Levison & Stern 2001; 
Bernstein et al. 2003). Theoretical models for the formation of KBOs have not addressed this 
issue. Theoretical predictions for r^, qs, and the space density of KBOs are also uncertain. 
Observations indicate that the observed space density of KBOs is / < 1% of the initial density 
of solid material in the planetesimal disk. Theoretical estimates of long-term collisional 
evolution yield / < 10% (e.g.. Stern & Colwell 1997b; Kenyon 2002), but the sensitivity of 
this estimate to the bulk properties of KBOs has not been explored in detail. 
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Here, we consider collision models for the formation and long-term collisional evolution 
of KBOs in the outer solar system. We use an analytic model to show how the break radius 
depends on the bulk properties and orbital parameters of KBOs (see Pan & Sari 2004, 
for a similar analytic model) and confirm these estimates with numerical calculations. If 
KBOs have relatively small tensile strengths and formed in a relatively massive solar nebula, 
we derive a break radius, r^ ~ 3-30 km, close to the observational limits. The numerical 
calculations also provide direct comparisons with the observed sky surface density and the 
total mass in the Kuiper Belt. Models with relatively weak KBOs and additional stirring 
by Neptune yield the best agreement with observations and make testable predictions for 
the surface density of KBOs at i? ~ 28-32. In the next 3-5 yr, occultation observations can 
plausibly test these predictions. 

We develop the analytic model in §2, describe numerical simulations of KBO evolution 
in §3, and conclude with a brief discussion and summary in §4. 



2. ANALYTIC MODEL 

2.1. Derivation 

We begin with an analytic collision model for the long-term evolution of an ensemble 
of KBOs. We assume that KBOs lie in an annulus of width Aa centered at a heliocentric 
distance a. KBOs with radius r orbit the Sun with eccentricity e and inclination i. The 
scale height H of the KBO population is H = a sin i. 

KBOs evolve through collisions and long range gravitational interactions. Here, we 
assume that gravitational interactions have reached a steady-state, with e and i constant in 
time. Based on numerical simulations of KBO formation, we adopt a broken power law for 
the initial size distribution, 

{^^j--as r < 1 km 

(1) 
fij^f-o'L r > 1 km 

with as = 3.5 for small objects and ai = 4.0 for large objects (Stern & Colwell 1997a; Davis 
& Farinella 1997; Stern & Colwell 1997b; Kenyon & Luu 1999a; Kenyon 2002; Kenyon & 
Bromley 2004). If V is the relative velocity, the collision rate for a KBO with radius ri and 
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all KBOs with radius r2 is 



'"^ ^ "^ ^ ^1 + 2.7^^) Vin + r,)\ (2) 



rft \4H a AaJ \ ' V^ 

where Ve is the escape velocity of a single body with mass, m = mi + m2 (Wetherill & 
Stewart 1993; Kenyon & Luu 1998, and references therein). 

The collision outcome depends on the impact velocity Vj and the disruption energy 
Qd- We follow previous investigators and define the energy needed to remove 50% of the 
combined mass of two colliding planetesimals, 

Qd = Qbr"" + pQgr"^ , (3) 

where Qtr^'' is the bulk (tensile) component of the binding energy and pQgV^^ is the gravity 
component of the binding energy (see, for example, Davis et al. 1985; Wetherill & Stewart 
1993; Holsapple 1994; Benz & Asphaug 1999; Housen & Holsapple 1999). The gravity com- 
ponent of the disruption energy varies linearly with the mass density p of the planetesimals. 
This expression ignores a weak relation between Qb,Qg and the impact velocity Vj 

g.,g.oc(^^^\ (4) 

where (3^ ~ 0.25-0.50 for rocky material and /3„ ^ —0.25 to —0.50 for icy material (e.g. 
Housen & Holsapple 1990, 1999; Benz & Asphaug 1999). For an analytic model with Vj = 
constant, the variation oi Qd with Vj is not important. We consider non-zero jS^ in complete 
evolutionary calculations described below. 



We adopt the standard center of mass collision energy (Wetherill & Stewart 1993) 

mim2V^ 

4(mi + 7712 



Q- = . """''j^.. ■ (5) 



where the impact velocity is 

Vf = V^ + Vl . (6) 

The mass ejected in a collision is 

vTiej = 0.5(mi + 1712) ( 77- ) , (7) 



Q 



where jSe is a constant of order unity (see Davis et al. 1985; Wetherill & Stewart 1993; 
Kenyon & Luu 1998; Benz & Asphaug 1999, and references therein). 
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For a single KBO, the amount of mass accreted in collisions with all other KBOs during 
a time interval 5t is 

/I 
— — 1112 dm2 . (8) 

The amount of mass lost is 

/I 
— — ruej dm2 . (9) 

KBOs with rhi > rha lose mass and reach zero mass on a removal timescale 

,,(,) « J!!£(lL . (10) 

mi -rua 

With these definitions, the evolution of an ensemble of KBOs depends on the relative 
velocity V , the size distribution, and the disruption energy. Because the disruption energy 
scales with size, larger objects are harder to disrupt than smaller objects. To produce a 
break in the size distribution, we need a 'break radius,' r?,, where 

tr < to for 1^ ^ f^b 

(11) 

tr > to for r > ri, 

and to is some reference time. We choose to = 1 Gyr as a reasonable e-folding time for the 
decline in the KBO space density. 



2.2. Application 

To apply the analytic model to the KBO size distribution, we use parameters appropriate 
for the outer solar system. We adopt i = e/2, with e = 0.04 for classical KBOs and e = 0.2 
for Plutinos. This simplification ignores the richness of KBO orbits, but gives representative 
results without extra parameters. We assume a total mass in KBOs, Mxbo, iii an annulus 
with Co = 40 AU and 6a = 10 AU. A minimum mass solar nebula has Mkbo ~ 10 M^; 
the current Kuiper Belt has Mkbo = 0.05-0.20 M® (Stern 1996a; Luu & Jewitt 2002; 
Bernstein et al. 2003). This range in initial KBO mass provides a representative range for 
the normalization constants, ns and n^, in our model size distribution. 

To model the destruction of KBOs, we adopt representative values for Qh, Qg, Ph, Pg, 
and p. Because the bulk properties of KBOs are poorly known, we consider wide ranges in 
Qb = 10^ to 10^ erg g -\ Cg = pQg = 10"^ to 10^ erg g-^, and f3g = 0.5-2.0. These ranges 
span analytic and numerical results in the literature (e.g., Davis et al. 1985; Holsapple 1994; 
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Love & Ahrens 1996; Benz & Asphaug 1999; Housen & Holsapple 1999). For simplicity, 
we adopt /Sb = 0; other choices have httle impact on the results. 

Figure 1 illustrates several choices for the disruption energy. We adopt a normalization 
constant for the gravity component of the disruption energy, 

Cg = pQ, = Co p ilO'y-''-^^ , (12) 

with p = 1.5 and Co = 1.5. All model curves then have the same gravity component of Qd 
at r = 1 km. The horizontal lines plot the impact energy for two massless KBOs with e = 
0.001 (lower line), 0.01 (middle line), and 0.1 (upper line). When e < 0.01, collisions cannot 
disrupt large KBOs with r > 1 km. Collisions disrupt smaller KBOs only if Qi, < 10® erg 
g^^ (e/0.01)^. When e ~ 0.1, collisions disrupt all small KBOs independent of their bulk 
strength. For KBOs with r > 1 km, disruptions are sensitive to /Sg and Cg. In general, more 
compact KBOs with larger f3g and Cg are harder to disrupt than fluffy KBOs with small Cg. 

Figure 2 plots several realizations of the analytic model as a function of KBO radius for 
a model with a mass equal to the minimum mass solar nebula and fragmentation parameters 
listed in the legend. For any combination of Qb and j3g, the removal timescale tr is a strong 
function of the KBO radius. Small KBOs with r < 1 km have large rhi/mQ and small removal 
timescales of 1 Myr or less. Larger objects lose a smaller fraction of their initial mass per 
unit time and have tr ~ 1-1000 Gyr. KBOs with smaller f3g are more easily disrupted and 
have smaller disruption times than KBOs with larger f3g. 

Figure 2 also illustrates the derivation of the break radius, r^,. The horizontal line at t = 
1 Gyr intersects the model curves at r^ ~ 3 km {Pg = 2) and at r^ ~ 6 km {Pg = 1.25). The 
break radius is clearly independent of the bulk strength Qb, and the reference time to, but is 
sensitive to f3g and e. At fixed collision energy, KBOs with smaller (3g fragment more easily. 
Thus, the break radius becomes larger as jSg becomes smaller. At fixed strength, KBOs with 
larger impact velocities also fragment more easily. Thus, the break radius becomes larger 
with larger e. 

As in Pan & Sari (2004), the break radius depends on the initial mass in KBOs and 
the collision velocity. For classical KBOs with e ~ 0.04, the break radius is r^ < 10-20 km 
when Mkbo ^ 10 M^ (Figure 3). In our model, Vb grows with the initial mass in KBOs. 
We derive roughly an order of magnitude change in r^ for a two order of magnitude change 
in the initial mass in KBOs. 

The break radius is also sensitive to Cg and (3g. KBOs with small Cg and (3g are easier to 
fragment than KBOs with large Cg and j3g. For a fixed mass in KBOs, plausible variations in 
Cg and j3g yield order of magnitude variations in r^ (Figure 3). At large Qb, these differences 
disappear. For a Kuiper Belt with 1% of the mass in a minimum mass solar nebula, the 
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model predicts no dependence of r?, on the physical variables for Qfc ^ 3 x 10^ erg g~^. This 
limit lies at Qi, > 10^ erg g~^ for a minimum mass solar nebula. 

Figure 4 shows how changes in the collision energy modify r^. Models with large eccen- 
tricity, e = 0.20, yield r^ ^ 3-100 km, compared to 0.3-10 km for e = 0.04. The predicted 
Tb is sensitive to the initial mass in KBOs, but is less sensitive to Cg and (3g. 

In contrast to Pan & Sari (2004), our analytic model indicates that a large break 
radius does not require a small bulk strength for the planetesimals. When the initial mass 
in KBOs is small, ~ 1% to 10% of the minimum mass solar nebula, the collision frequency 
for 10-100 km KBOs is also small, ~ 10-100 Gyr^^. To remove sufficient material in 5 Gyr, 
these collisions must produce mostly debris. Thus, Qd must be small. Large bulk strengths, 
Qb ^ 10^ erg g~^, preclude debris-producing collisions (Figure 1) and result in small Vf, 
(Figures 3-4). As the initial mass in KBOs increases to 10% to 100% of the minimum 
mass solar nebula, collision frequencies also increase. More frequent collisions can remove 
large KBOs from the size distribution even when Qt is large. Thus, the break radius is less 
sensitive to Qb for massive nebulae. 



2.3. Implications for the Kuiper Belt 

The conclusions derived from the analytic model have direct implications for the for- 
mation and evolution of KBOs. The apparent break in the sky surface density at i? ~ 28 
requires a break in the size distribution at r?, ^ 20 km for an albedo of 0.04. The analytic 
model shows that a measured r^ > 20 km requires a nebula with a mass in solids of at least 
10% of the minimum mass solar nebula. Producing such a large break radius from collisions 
is easier in an initially more massive nebula. The current mass in KBOs is < 1% of the 
minimum mass solar nebula. Thus, the analytic model provides additional support for an 
initially more massive solar nebula at 30-50 AU (see also Stern & Colwell 1997a; Kenyon 
& Luu 1999a; Kenyon 2002, and references therein). 

The break in the sky surface density also favors a low bulk strength, Qb < 10^ erg g~^, 
for KBOs. This Qb is smaller than the Qb ~ 10^ erg g~^ derived from numerical models 
of collisions of icy bodies (e.g., Benz & Asphaug 1999). However, a low bulk strength is 
consistent with the need for a strengthless rubble pile in models of the break-up of comet 
Shoemaker-Levy 9 (e.g., Asphaug & Benz 1996). 

The analytic model may also explain differences in the observed size distributions for 
different dynamical classes of KBOs. From Figures 3-4, resonant KBOs and scattered KBOs 
with large e and i should have a larger rf, than classical KBOs with small e but large i (Luu & 



Jewitt 2002). The group of 'cold' classical KBOs with small e and small i (Levison & Stern 
2001) should have even smaller r^. The observations provide some support for this division. 
Relative to the number of bright KBOs, there are fewer faint resonant and scattered KBOs 
and more classical KBOs than expected (Bernstein et al. 2003). Long-term collisional 
evolution could be responsible for removing higher velocity resonant and scattered KBOs 
and leaving behind lower velocity classical KBOs. 

To test the analytic model and provide direct comparisons with the observations, we 
need a numerical model of accretion and erosion. A numerical model accurately accounts 
for the time dependence of the total mass in KBOs and thus provides a clear measure of the 
removal time for a range of sizes. Numerical models also yield a direct calculation of the size 
distribution and thus measure the 'depletion' of KBOs as a function of radius. 



3. NUMERICAL MODELS 

To test the analytic model, we examine numerical simulations with a multiannulus 
coagulation code (Kenyon & Bromley 2004, and references therein). For a set of A^ concentric 
annuli surrounding a star of mass M, this code solves numerically the coagulation and 
Fokker-Planck equations for bodies undergoing inelastic collisions, drag forces, and long- 
range gravitational interactions (Kenyon & Bromley 2002). We adopt collision rates from 
kinetic theory and use an energy-scaling algorithm to assign collision outcomes (Davis et al. 
1985; Wetherill & Stewart 1993; Weidenschilling et al. 1997; Benz & Asphaug 1999). We 
derive changes in orbital parameters from gas drag, dynamical friction, and viscous stirring 
(Adachi et al. 1976; Ohtsuki, Stewart, & Ida 2002). The appendix describes updates to 
algorithms described in Kenyon & Bromley (2001, 2002a, 2004) and Kenyon & Luu (1998, 
1999). 

We consider two sets of calculations. For models without gravitational stirring, we set 
e = constant, i = e/2, and calculate the collisional evolution of an initial power law size 
distribution. These calculations require a relatively small amount of computer time and 
allow a simple test of the analytic model. 

Complete evolutionary calculations with collisional evolution and gravitational stirring 
test whether particular outcomes are physically realizable. These models also provide direct 
tests with observables, such as the current mass in KBOs and the complete KBO size distri- 
bution. Because these models do not allow arbitrary e and i, they are less flexible than the 
constant e models. 

To provide some flexibility in models with gravitational stirring, we calculated models 
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with and without stirring by Neptune at 30 AU. In models without Neptune, large KBOs 
with radii of 1000-3000 km stir up smaller KBOs to the disruption velocity. The KBO size 
distribution, including the break radius, then depends on the radius of the largest KBO 
fL,KBO formed during the calculation. Because tlkbo depends on Qi, and Qg (e.g., Kenyon 
&; Luu 1999a; Kenyon 2002), r^ also depends on Qb and Qg. 

In models with Neptune, long-range stirring by Neptune can dominate stirring by local 
large KBOs. The break radius then depends on the long-range stirring formula (Weiden- 
schilling 1989; Ohtsuki, Stewart, & Ida 2002) and the timescale for Neptune formation. 
Here, we assume a 100 Myr formation time for Neptune, whose semimajor axis is fixed at 
30 AU throughout the calculation. The mass of Neptune grows with time as 

' 6 X 10^7 e(*-*o)/*i g t < to 

MMep ^ <! 6 X 10^7 g + CNep{t - h) to < t < ta (13) 

1.0335 X 10^9 g t>t2 

where CiVep is a constant and toj ^i) and t2 are reference times. For most calculations, we set 
to = 50 Myr, ti = 3 Myr, and t2 = 100 Myr. These choices allow our model Neptune to reach 
1 M0 in 50 Myr, when the largest KBOs have formed at 40-50 AU, and reach its current 
mass in 100 Myr. This prescription is not intended as a model for Neptune formation, but 
it provides sufficient extra stirring to test the prediction that the break radius depends on 
the amount of local stirring. 



3.1. Constant eccentricity models 

Calculations with constant eccentricity allow a direct test of the analytic model. We 
performed a suite of ~ 200 4.5 Gyr calculations for a range in fragmentation parameters, 
with log Qh = 1-8, Pg = 0.5-2.0, and log Cg + 5 Pg = 0.01-20. The initial size distribution 
of icy planetesimals has sizes of 1 m to 100 km in mass bins with 6 = mi^i/mi = 1.7 and 
equal mass per mass bin. The planetesimals lie in 32 annuli extending from 40 AU to 75 AU. 
The central star has a mass of 1 Mq. The initial surface density. So = 10^'^ g cm^^ to 10^^ 
g cm~^ at 40 AU, ranges from 1% to 100% of the minimum mass solar nebula extended to 
the Kuiper Belt (Weidenschilling 1977; Hayashi 1981). The initial eccentricity, e = 0.04 
and 0.20, spans the observed range for classical and resonant KBOs (Luu & Jewitt 2002). 

In Kuiper Belt models with large e, fragmentation is the dominant physical process (see 
also Stern & Colwell 1997a; Kenyon & Luu 1999a; Kenyon & Bromley 2002). Large, 10- 
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100 km objects grow very slowly. Smaller objects suffer numerous disruptive collisions that 
produce copious amounts of debris. Debris fills lower mass bins, which suffer more disruptive 
collisions. In 10-100 Myr, this collisional cascade reduces the population of 0.1-1 km and 
smaller objects. The size distribution then follows a broken power law, with a^ ^ 3 for large 
objects and 0:5 ~ 2.5 for the small objects (Dohnanyi 1969; Williams & Wetherill 1994; 
Pan & Sari 2004). 

As the evolution proceeds, the size distribution evolves into a standard shape (Figure 
5). After 4.5 Gyr, the largest objects grow from 100 km to 125 km. From ~ 10 km to 125 
km, the size distribution continues to follow a power-law with a^ ~ 3. Test calculations 
suggest that this power law slope is fairly independent of the initial power law. 

At smaller sizes, the shape of the size distribution depends on the bulk strength. For 
Qb ^ 10^ erg g~^, disruptive collisions produce a break in the size distribution; the power 
law slope changes from aL ^ S to as ^ 2.5. For Qi, < 10*^ erg g~^, disruptive collisions 
produce two breaks, one at 1-10 km and another at ~ 0.1 km (Figure 5). The break at 1-10 
km is where growth by accretion roughly balances loss by disruption. The break at ~ 0.1 km 
is where debris produced by collisions of larger objects roughly balances loss by disruptive 
collisions. Between these two sizes, the slope of the size distribution ranges from aj ~ —0.5 
for Qb ~ 10 erg g~^ to «/ ~ 1 for Qb ~ 10^ erg g~^. The slope of the power law for small 
sizes ranges from aj ~ 4.5-5.0 for Qi, ~ 10 erg g~^ to «/ ^ 3 for Qb ~ 10^ erg g~^. As Qb 
approaches 10^ erg g~^, the slopes of both power laws converge to a ~ 2.5. 

We define the first inflection point in the size distribution as the break radius. To 
measure r^, we use a least-squares fit to derive the best-fitting power-law slopes, a^ and 
ai, to the calculated size distribution. Using Poisson statistics to estimate errors in A^, we 
derive r^ and its la error by minimizing the residuals in the fits. Typical errors in log Vb are 
± 0.02-0.05. Tests indicate that the derived r^ is more sensitive to the mass resolution of our 
calculations, 6, than to the range in log r used in the fits. This error is also small compared 
to the range in log Vb, ~ 0.2, derived from repeat calculations with the same combination of 
Qb and jSg. 

Figure 6 shows results for calculations with constant e = 0.04. For all initial masses 
and Qb < 10^, r^ is independent of Qb- For small initial masses, calculations with stronger 
objects yield small break radii. At larger initial masses, this sensitivity to Qb disappears. 
Calculations for minimum mass solar nebulae show little variation of Vb with Qb- 

These results confirm the basic features of the analytic model. Both models predict Vb < 
1 km for low mass nebulae with ~ 1% of the mass in the minimum mass solar nebula. Larger 
break radii, ~ 1-10 km, are possible in more massive nebulae. The analytic model predicts 
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large break radii for larger e than the numerical calculations. In numerical calculations with 
e = 0.2, disruptive collisions reduce the space density considerably in ~ 100 Myr. The 
smaller collision rates prevent formation of a break in the size distribution at large radii. 
Thus, numerical calculations with e = 0.2 yield log r^ only ~ 0.1-0.2 larger than calculations 
with e = 0.04. 



3.2. Full evolution models 

These calculations begin with 1-1000 m planetesimals in mass bins with 5 = 1.4 or 
1.7 and equal mass per bin. The planetesimals lie in 32 annuli at 40-75 AU. Models with 
Neptune have an extra annulus at 30 AU. For most models, we adopt Cq = 10"^ or cq = 10~^ 
and i = e/2 for all planetesimals. At the start of our calculations, these initial values yield a 
rough balance between viscous stirring by 0.1-1 km objects and collisional damping of 10- 
100 m objects. The bodies have a mass density pd = 1.5 g cm~^, which is fixed throughout 
the evolution. We consider a range in initial surface density, with Sq = 0.03-0.3 g cm~^ 
(ao/30 AU)-3/2. 

To measure the sensitivity of our results to stochastic variations, we performed 2-5 
calculations for each set of fragmentation parameters. For /5fe = and a factor of ten range 
in So, we considered log Qi, = 1, 2, 3, 4, 5, 6, and 7; Cq = 0.15 and 1.5; and j3g = 1.25 and 
2.0. We also performed a limited set of calculations for Cg = 0.15 and 1.5, f3g = 0.5, and a 
small set of log Qb- Although stochastic variations can change the size of the largest object 
at 40-50 AU, repeat calculations with identical initial conditions yield small changes in the 
shape of the size distribution or the location of the break radius. A few calculations with 
Pb ¥" ^ yield interesting behavior in the size distribution at 1-100 m sizes, but r^ does not 
change dramatically. A larger suite of calculations with /5„ 7^ leads to similar conclusions. 
We plan to report on these aspects of the calculations in a separate paper. 

Icy planet formation in the outer solar system follows a standard pattern (see Kenyon 
& Luu 1999a; Kenyon & Bromley 2004; Goldreich, Lithwick, & Sari 2004). Small plan- 
etesimals with Ti < 1 km first grow slowly. CoUisonal damping brakes the smallest objects. 
Dynamical friction brakes the largest objects and stirs up the smallest objects. Gravitational 
focusing factors increase, and runaway growth begins. At 40-50 AU, it takes ~ 1 Myr to 
produce 10 km objects and another 3-5 Myr to produce 100 km objects. Continued stirring 
reduces gravitational focusing factors. Collisions between the leftover planetesimals produce 
debris instead of mergers. Runaway growth ends and the collisional cascade begins. 

During the collisional cascade, the mass in 1-10 km and smaller objects declines precipi- 
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tously. Because gravitational focusing factors are small, collisions between two planetesimals 
are more likely than collisions between a planetesimal and a 100-1000 km protoplanet. Thus, 
disruptive collisions grind leftover planetesimals into small dust grains, which are removed by 
radiation pressure and Poynting-Robertson drag (e.g.. Burns, Lamy, & Soter 1979; Takeuchi 
& Artymowicz 2001). At 40-50 AU, the surface density falls by a factor of two in 100-200 
Myr, a factor of 4-5 in 1 Gyr, and more than an order of magnitude in 3-4 Gyr (see also 
Kenyon & Bromley 2004). After 4.5 Gyr, the typical amount of solid material remaining at 
40-50 AU is 3% to 10% of the initial mass (see below). 

As collisions and radiation remove material from the system, the largest objects continue 
to grow slowly. In most calculations, it takes 10-50 Myr to produce the first 1000 km object. 
The largest objects then double their mass every 100 Myr to 1 Gyr. After 4.5 Gyr, the 
largest objects have radii ranging from ~ 100 km {Qi, < 10^ erg g~^, f3g = 0.5) to 5000 km 
{Qb > 10^ erg g-\ I3g = 2.0). Calculations with f3g = 1.25 and Qb ~ 10^ erg g-^ to 10^ 
erg g^^ favor the production of objects with radii of 1000-2000 km, as observed in the outer 
solar system. 

These general results are remarkably independent of the initial conditions and of some 
input parameters (Kenyon & Luu 1999a; Kenyon & Bromley 2004). The size distribution 
of objects remaining at 4.5 Gyr is not sensitive to the initial disk mass, the initial size 
distribution, the initial eccentricity and inclination (for cq ^ 10~^), the mass resolution 6, 
the width of an annulus 6a, or the gas drag parameters. The orbital period P and the surface 
density set the collision timescale, t oc P/S. Although stochastic variations can produce 
factor of 2 or smaller variations in growth times, all timescales depend on the collision time 
and scale with the current surface density. 

To derive r^ for these calculations, we again examine the size distribution at 4.5 Gyr 
(Figure 7). For radii Tj ~ 10-1000 km, the size distribution follows a power-law, N oc r~"^ 
with slope aL ~ 3-3.5. For smaller radii, the size distribution has inflection points at log 
Tj ^ — 1 to 1 and at log r^ < —1. Between the two inflection points, the power-law slope is 
shallow, with aj ~ 0-2. For small objects, the power-law slope is steep, with 0:5 ^ 2-4. 

The slope of the intermediate power-law depends on the fragmentation parameters. 
Calculations with small Qb and Pg yield small a/. For larger Qb and Pg, the slope approaches 
the collisional limit, aj ~ 2.5 (Dohnanyi 1969; Williams & Wetherill 1994). The extent 
of the shallow, intermediate power-law also varies with Qb and f3g. For Qb < 10'^ erg g~^ 
and f3g ~ 1.25, the size distribution is relatively fiat from log Vi ^ —1 to 0.0-0.5 (Figure 7). 
Large Qb and /Sg result in a smaller extent for the intermediate power law. 

Figure 8 shows results for r^ as a function of log Qb for several sets of calculations. All 
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calculations began with a mass in solids comparable to a minimum mass solar nebula. For 
models without stirring by Neptune, we derive log rf, ^ 5.0-5.6 (log Qi, < 4) and log rf, ^ 
5. 3 ±0.3 — 0.5 log Qb (log Qb > 4). Although there is some overlap in the results for different 
sets of parameters, a weaker gravity component to the disruption energy {(3g ^ 1.25) favors 
larger r^. This result confirms the conclusion derived from the analytic model. 

Stirring by Neptune yields larger values for the break radius (Figure 8). As Neptune 
reaches its final mass at 80-100 Myr, long-range stirring rapidly increases the eccentricities 
of objects at 40-50 AU. This stirring accelerates the collisional cascade, which depletes the 
population of small planetesimals and halts the growth of the largest objects. The larger 
e and longer duration of the collisional cascade moves the break in the size distribution to 
larger radii. 

Figure 9 shows size distributions at 4.5 Gyr for three models with stirring by Neptune. 
The legend lists input values for Qg and log Qb- In these calculations, the break is at rj, ^ 
5-10 km, compared to r^ ^ 1-5 km for models without Neptune stirring. The position of 
the break is relatively independent of Qb or Qg (see Figure 8). For KBOs with sizes smaller 
than the break, the slope of the size distribution is sensitive to Qb but not to Qg. Models 
with Qb > 10^ erg g~^ have more small objects with radii of ~ 0.1 km than models with 
g,<10^ergg-i. 



3.3. Comparisons with observations of KBOs 

The analytic model and the numerical evolution calculations yield a consistent picture 
for the size distribution of icy bodies in the outer solar system. The general shape of the size 
distribution does not depend on the initial conditions or input parameters. The typical size 
distribution has two power laws - one power-law for large objects with radii > 10-100 km 
and a second power-law for small objects with radii < 0.1 km - connected by a transition 
region where the number of objects per logarithmic mass bin is roughly constant. The power 
law slope for the large objects is also remarkably independent of input parameters and initial 
conditions. 

The fragmentation parameters and the amount of stirring set the location of the tran- 
sition region and the power-law slope for the small objects. In our calculations, the break 
radius is Vb > 10 km when icy objects are easy to break and the stirring is large. Strong 
icy objects and small stirring favor a small break radius, t^ < 1 km. When the break radius 
is small, the extent of the transition region is also small, less than an order of magnitude 
in radius. When the break radius is large, the transition region can extend for 2 orders of 
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magnitude in radius. 

To compare the numerical results to observations of KBOs, we convert a calculated 
size distribution into a predicted sky surface density of KBOs as a function of apparent 
magnitude. Current observations suggest that the observed sky surface density, a{m) - the 
number of KBOs per square degree on the sky, follows a power law 

log (T(m) = log (To + a{m — mo) (14) 

where a, ao and mo are constants. With a ~ 0.6 and mo = 23, this function fits the 
data fairly well for R-band magnitudes, R ~ 22-26. For i? < 22 and i? > 26, the simple 
function predicts too many KBOs compared to observations (Bernstein et al. 2003). The 
observations also suggest that the scattered and Plutino populations of KBOs have different 
surface density distributions than classical KBOs, with a ~ 0.6 for scattered KBOs and 
Plutinos and a ~ 0.8 for classical KBOs (Bernstein et al. 2003). 

Because we do not include the dynamics of individual objects in our calculations, we 
cannot predict the relative numbers of KBOs in different dynamical classes. However, we 
can predict a{m) for all KBOs and see whether the calculations can explain trends in the 
observations. 

To derive a model cr(m), we assign distances (Iq to an ensemble of objects chosen ran- 
domly from the model size distribution. For a random phase angle (3 between the line-of-sight 
from the Earth to the object and the line-of-sight from the Sun to the object, the distance 
of the object from the Earth is dE = (IqCosP — (1 + (i|(cos^/3 — 1))^^^. The red magnitude 
of this object is i? = -Rq + 2.5 log (ti/t2) — 5 log rxBO, where Ro is the zero point of the 
magnitude scale, r^ is the radius of the object, ti = 2dQdE, and t2 = ^((1 — g)'Pi + 5*02) 
(Bowell et al. 1989). In this last expression, u is the albedo, and g is the slope parameter; 
01 and 02 are phase functions that describe the visibility of the illuminated hemisphere of 
the object as a function oi f3. We adopt standard values, uj = 0.04 and g = 0.15, appropriate 
for comet nuclei (Jewitt et al. 1998; Luu & Jewitt 2002; Brown & Trujillo 2004). The 
zero point Rq is the apparent red magnitude of the Sun, mpt^Q = —27.11, with a correction 
for the V-R color of a KBO, Ro = ttiji^q + 5(V-R)i^Bo- Observations suggest that KBOs 
have colors that range from roughly —0.1 to 0.3 mag redder than the Sun (Jewitt & Luu 
2001; Tegler & Romanishin 2003; Tegler, Romanishin, & Consolmagno 2003). We treat 
this observation by allowing the color to vary randomly in this range. 

For this application, we assign distances c/q = 40-50 AU and derive the number of 
objects in half magnitude bins for R = 15-50. In most models, the surface density of objects 
predicted by the model closely follows the linear relation for log a with a ~ 0.55-0.7 and 
Ro ~ 21-24 (see also Kenyon & Luu 1999b; Kenyon 2002). To make easier comparisons 
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between observations and theory, we follow Bernstein et al. (2003) and define the relative 
space density as the ratio between the model and the linear surface density relation with a 
= 0.6 and Rq = 23. 

Figure 10 shows the relative surface density for several KBO calculations. For bright 
KBOs with R ^ 22-28, the surface density closely follows the linear relation, equation (14). 
This result is independent of the fragmentation parameters, the initial mass at 40-50 AU, 
the initial size distribution, and the amount of stirring by Neptune. Thus, the coagulation 
models provide a robust prediction for a at i? ~ 22-28 (see also Kenyon & Luu 1999b; 
Kenyon 2002). 

There are significant differences between the models at i? < 22 and R > 28. For R < 
22, the models fall into two broad classes defined by the ratio of the disruption energy to 
the typical collision energy, Qd/Qi- Calculations with large Qd/Qi produce many large, 
bright KBOs. Calculations with small Qd/Qi produce few large, bright KBOs. The range 
in production rates is roughly 4 orders of magnitude. 

For fainter KBOs, the differences between models become even more significant. At R ^ 
27-30, collisions remove weak KBOs from the size distribution. Thus, most calculations with 
Neptune stirring exhibit a drop in the relative surface density. Calculations with Qh < 10^ 
erg g~\ Qg < 0.1-0.2, and no Neptune stirring also produce fewer KBOs at i? ~ 27-30. 

Relative to the nominal power-law, all calculations produce an excess of KBOs at i? ~ 
29-33. For weak KBOs with Qi, < 10^ erg g~^, the model predicts a factor of 3 excess 
compared to the power-law. For strong KBOs with Qi, > 10^ erg g"\ this excess grows to a 
factor of 10-30. Stirring by Neptune has little impact on this excess. 

For fainter KBOs with R > 32, stirring by Neptune is very important. Our models 
produce a 5-12 order of magnitude deficit of KBOs relative to the power-law surface density 
relation. Weaker KBOs produce larger deficits. Neptune stirring also produces larger deficits. 
We derive the largest deficit, 12 orders of magnitude at i? ~ 40, for models with Neptune 
stirring, Q^ < 10'^ erg g~^, and Qg < 0.15 erg g^^ (Figure 10). 

The derived a{m) yields good agreement with the observations (Figure 11). For R ^ 
22-28, most models account for the variation of relative number density with R magnitude. 
Calculations with weaker KBOs, Qb ^ 10^ erg g~^, reproduce the dip in the relative number 
density at /? ~ 19-20. Our ensemble of models suggests that the magnitude of the dip 
depends more on stochastic phenomena than on model parameters. The small Q^ models 
also provide better agreement with observations for fainter KBOs with R ^ 26-28. At 
R ^ 20-21 and at i? ~ 29-30, models without Neptune stirring produce an excess of KBOs 
relative to the observations; models with Neptune stirring yield better agreement with the 
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data. Thus, models with weak KBOs, Qb ^ 10^ erg g ^, and with Neptune stirring provide 
the best explanation for current observations of the shape of the size distribution of KBOs. 

In addition to the relative size distribution, the collision models provide fair agreement 
with the absolute numbers of KBOs (Figure 12). Current data suggest a total mass of ~ 
0.1 Me in KBOs at 35-50 AU (Luu & Jewitt 2002; Bernstein et al. 2003). To form KBOs 
by coagulation in 10-100 Myr, collision models require an initial mass in solids comparable 
to the minimum mass solar nebula. This result suggests that the current mass in KBOs is 
~ 1% of the initial mass in solid material at 35-50 AU. After 4.5 Gyr, our collision models 
have 3%-10% of the initial mass in 1 km and larger objects. Models with Neptune stirring 
are more efficient at removing material from the size distribution (Figure 12). 

This result is encouraging. Once large objects form in the Kuiper Belt, the collisional 
cascade can remove almost all of the leftover planetesimals, which contain 90% to 97% of 
the initial mass (see also Stern & Colwell 1997a; Kenyon & Luu 1999a). Other processes 
not included in our calculations will also remove large objects. Dynamical interactions with 
Neptune and other giant planets can remove 50% to 80% of the initial mass (Holman & 
Wisdom 1993; Levison & Stern 1995). Interactions with field stars can also remove KBOs 
(Ida et al. 2000). Including these processes in a more realistic collision calculation should 
bring the predicted number of KBOs into better agreement with observations. We plan to 
describe calculations testing the role of Neptune in the formation and evolution of the Kuiper 
Beh. 

Although we cannot develop tests using data for KBO dynamical families, the models 
provide some insight into general trends of these observations. Because scattered KBOs 
and most resonant KBOs have had close dynamical interactions with Neptune, these objects 
probably formed closer to the Sun than classical KBOs. If Neptune stirring halted accretion 
in the Kuiper Belt, this difference in heliocentric distance can produce an observable differ- 
ence in KBO sizes. For a formation timescale, t oc P/S oc a^, KBOs at 45 AU take twice as 
long to form as KBOs at 35 AU. During the late stages of runaway growth, this difference 
in formation timescales leads to a factor of ~ 2 difference in the maximum size of a KBO. 
Although the oligarchic growth phase erases this difference, significant stirring by Neptune 
during runaway growth might preserve the difference and lead to the apparent lack of large 
classical KBOs (formed at ~ 45 AU) relative to resonant KBOs (formed at ~ 35 AU). We 
plan additional numerical calculations to test this idea. 
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4. CONCLUSIONS AND SUMMARY 

We have developed an analytic model for the formation of a break in the power law 
size distribution of KBOs in the outer solar system. For a mass in KBOs equivalent to the 
current mass at 40-50 AU and for e = 0.04 orbits, the model predicts a break at a radius, 
Th ~ 0.1-1 km. For a massive Kuiper Belt with e = 0.2, the break moves to r^ ~ 10-100 km. 
These results agree with the model of Pan & Sari (2004). 

In contrast to Pan & Sari (2004), our model predicts a smaller sensitivity to the bulk 
strength of KBOs. When the mass in KBOs is 1% of the minimum mass solar nebula, r^ is 
independent of Qb for Qb < 10^ erg g^^. For Qi, > 10^ erg g^^, r?, declines with increasing 
Qb- As the total mass in KBOs increases, the break radius is less sensitive to Qb- When the 
mass in KBOs is comparable to a minimum mass solar nebula, the analytic model suggests 
that Tb is roughly constant for all reasonable Qb- 

To test the analytic model, we used a suite of numerical simulations. Constant eccen- 
tricity calculations, with no velocity evolution due to gravitational interactions, confirm the 
analytic results. For models with constant e = 0.04-0.20, the break radius depends on Qb 
and the total mass in KBOs, 

1-0 iokY^" {w^y^ km Qb < 10^ erg g"^ 

Tbr- ~ < (15) 



^■Hokf"{w^) km g.>10^ergg-i 



where Mkbo,o is the mass of a minimum mass solar nebula extended into the Kuiper Belt, 
~ 10 Me at 40-50 AU. 

Simulations of complete KBO evolution, with velocity stirring, generally require more 
initial mass in planetesimals to yield the same results. In models without Neptune formation, 
stirring by large objects with radii of 1000-3000 km yield 

f 1 - 3 km Qb< 10^ erg g-^ 

Tbr^ I 1/2 ^^^^ 

[ 1 - 3 (tot^-t) km Qb > 10^ erg g-^ 

for models starting with a mass in solid comparable to the minimum mass solar nebula. 
Calculations with Neptune at 30 AU allow larger break radii independent of Qb, with 

rbr^ 3 -10 km g5<10^ergg-i (17) 

In both cases, models with more initial mass yield larger r^. 
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Comparisons between observed and predicted size distributions of KBOs allow tests of 
models for KBO formation and evolution. For a broad range of input parameters, KBO 
models with and without stirring by Neptune yield good agreement with observations for 
R ^ 21-27. The observed surface density of brighter KBOs suggests that KBOs have 
Qb ^ 10^ — 10^ erg g~^. Although stirring by Neptune modifies the shape of the KBO size 
distribution for i? < 22, the observations do not discriminate clearly between models with 
and without Neptune. Improved observational constraints on the surface density of KBOs 
for R < 21-22 might provide tests for the relative formation times of Neptune and large 
KBOs. 

Observations for i? > 27 may also yield constraints on the formation of Neptune. Long- 
range stirring by Neptune is more important for the size distribution of fainter KBOs, R > 
27. Most models predict a small dip in the size distribution at i? ~ 27-30, a peak at i? ~ 
30-34, and a deep trough at i? ~ 35-45. Because stirring by Neptune removes more objects 
with radii of 1-30 km from the KBO size distribution, models with Neptune produce a larger 
dip and a deeper trough than models without Neptune. The depth of the small dip in models 
with Neptune stirring is close to the depth observed in recent HST observations (Bernstein 
et al. 2003). 

Calculations with Neptune stirring yield many orders of magnitude fewer KBOs with 
R ^ 32-42 than calculations without Neptune. Current observations do not probe this 
magnitude range. However, ongoing and proposed campaigns to detect small KBOs from 
occultations of background stars allow tests of the models (Bailey 1976; Brown & Webster 
1997; Roques & Moncuquet 2000; Cooray & Farmer 2003; Roques et al. 2003). For u = 
0.04, detections of 1 km KBOs provide constraints at i? ~ 33-35, where Neptune stirring 
models predict a sharp drop in the KBO number density. Direct detections of smaller KBOs, 
with radii of ~ 0.1 km, constrain model predictions at i? ~ 40, where Neptune stirring models 
predict 3-6 order of magnitude fewer KBOs than models without Neptune stirring. 

After 4.5 Gyr of collisional evolution, all of the numerical calculations predict a small 
residual mass in large KBOs. For Q^ < 10® erg g"^, the simulations leave / ~ 3%-8% of the 
initial planetesimal mass in KBOs with radii of 1 km and larger. Models with stirring by 
Neptune contain less mass in KBOs than models without Neptune. Calculations with Qf, > 
Xq6 gj,g g-i iiave a larger range in /; models with stirring by Neptune still leave ~ 3%-5% 
of the initial mass in large KBOs. 

These results are encouraging. Although our calculations leave more material in large 
objects than the current mass in KBOs, ~ 0.5%-l% of a minimum mass solar nebula, other 
processes can reduce the KBO mass considerably. Formation of Neptune at 10-20 Myr, 
instead of our adopted 100 Myr, probably reduces our final mass estimates by a factor of 
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two. Dynamical interactions with Neptune and passing stars can also remove substantial 
amounts of material (e.g., Holman & Wisdom 1993; Levison & Stern 1995; Duncan et 
al. 1995; Malhotra 1996; Levison & Duncan 1997; Morbidelli & Valsecchi 1997; Ida et 
al. 2000). These studies suggest that a combination of collisional grinding and dynamical 
interactions with Neptune or a passing star can reduce a minimum mass solar nebula to 
the mass observed today in the Kuiper Belt. We plan to describe additional tests of these 
possibilities in future publications. 

Finally, our calculations provide additional evidence that observations of Kuiper Belt 
objects probe the formation and early evolution of Neptune and other icy planets in the 
outer solar system. Better limits on the sizes of the largest KBOs probe the timescale for 
Neptune formation^. These observations also constrain the bulk strength of KBOs during 
the formation epoch. The detection of small KBOs, r ^ 0.1-10 km, by occultations (e.g. 
TAOS; Marshall et al. 2003) or by direct imaging (e.g., OWL; Gilmozzi, Rierckx, & Mon- 
net 2001) yields complementary constraints. As the observations improve, the theoretical 
challenge is to combine collisional (this paper; Goldreich, Lithwick, & Sari (2004)) and 
dynamical (e.g. Malhotra 1995; Gomes 2003; Levison & Morbidelli 2003; Quillen, Trilling, 
& Blackman 2004) calculations to derive robust predictions for the formation and evolution 
of Uranus, Neptune, and smaller icy planets at heliocentric distances > 15 AU. Together, 
the calculations and the observations promise detailed tests of theories of planet formation. 



We acknowledge a generous allotment, ~ 3000 cpu days from February 2003 through 
March 2004, of computer time at the supercomputing center at the Jet Propulsion Laboratory 
through funding from the NASA Offices of Mission to Planet Earth, Aeronautics, and Space 
Science. Advice and comments from M. Geller and S. A. Stern improved our presentation. 
We thank P. Michel and A. Morbidelli for extensive discussions that improved our treatment 
of collisional disruption. We also acknowledge discussions with P. Goldreich and M. Holman. 
The NASA Astrophysics Theory Program supported part of this project through grant NAG5- 
13278. 



A. APPENDIX 

Kenyon & Luu (1998, 1999) and Kenyon & Bromley (2001, 2002a, 2004) describe 



^The discovery of Sedna (Brown, Trujillo, & Rabinowitz 2004) tests models for KBO formation at 50-100 
AU. 
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algorithms and tests of our multiannulus planet formation code. Here, we describe an update 
to the fragmentation algorithm. 

In previous calculations, we used fragmentation prescriptions summarized by Davis et 
al. (1985) and Wetherill & Stewart (1993). Both methods write the strength S* of a pair of 
colliding objects as the sum of a constant bulk strength, 5*0, and the gravitational binding 
energy. Eg, 

S = So + Eg. (Al) 

For Eg oc {rrii + mj)/rij, the strength is S* = 5*0 + Sirfj, where Si is a constant and Vij is the 
radius of an object with mass rrii + mj. When the collision energy, Q/, exceeds S, Wetherill 
& Stewart (1993) derive the mass lost by catastrophic disruption as the ratio of the impact 
energy to a crushing energy Q^. Davis et al. (1985) assume that a fixed fraction, /x_b, of 
the impact kinetic energy is transferred to ejected material and derive the fraction of the 
combined mass lost to disruption. In most cases, the Davis et al. (1985) algorithm yields 
less debris than the Wetherill & Stewart (1993) algorithm. 

To take advantage of recent advances in numerical simulations of collisions (e.g., Benz 
& Asphaug 1999; Michel et al. 2001, 2002, 2003), we now define a disruption energy 
Qd required to eject 50% of the combined mass of two colliding bodies (equation 3 in the 
main text). For collision energy Q/, the mass ejected in a catastrophic collision is m^j = 
0.5(mi + m2){Qi/QdY^- For most applications we set (3e = 1.125 (e.g., Davis et al. 1985). 
When mej < 10~^ {rrii + rrij), we follow Wetherill & Stewart (1993) and set rriej = 0. 

To derive the size and velocity distribution of the ejected material, we adopt a simple 
procedure for all collisions. We define the remnant mass. 

The mass of the largest ejected body is 

mL,ej = 0.2mej (A3) 

We adopt a cumulative size distribution for the remaining ejected bodies, ndTn) oc m^^, with 
b = 0.8 (Dohnanyi 1969; Williams & Wetherill 1994), and require that the mass integrated 
over the size distribution equal TJiej- We assume that all bodies receive the same kinetic 
energy per unit mass, given by the initial relative velocities of the two bodies, 

VS = hl + vl + h] + v] , (A4) 

where h and v are the horizontal and vertical components of the velocity dispersion relative 
to a circular orbit (e.g. Kenyon & Luu 1998). We derive mass- weighted h!^ and v\ for the 
combined object and the debris, V^j = (h[Y + {v[Y. 
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This procedure, which we apply to cratering and disruptive colhsions, is computation- 
ally efficient and maintains the spirit of recent analytic models and numerical simulations. 
Comparisons with our previous results using the Davis et al. (1985) and Wetherill & Stewart 
(1993) algorithms suggest that the new algorithm yields intermediate 'mass loss rates' for Pg 
= 2 and Qb ~ 10^-10^ erg g^^ Calculations with /3g ~ 1.2-1.5 (Benz & Asphaug 1999) yield 
larger mass loss but do not change the results significantly. When the bulk strength depends 
on the particle radius, the size distribution for small objects with Tj < 0.1 km depends on 
the exponent f3b of the bulk component of the strength. We plan to report on the details of 
these differences in future papers on the formation of KBOs and terrestrial planets. 
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Fig. 1. — Comparison of catastrophic disruption energy for models described by equation 
(3) with p = 1.5 g cm~^. Sohd curve: Qb = 1.6 x 10'' erg g~^, Qg = 1.5 erg g~\ Pb = —0.42, 
f]g = 1.25 (Benz & Asphaug 1999). Dot-dashed curves: Qg = 2.67 x 10"^ erg g-^ f3b = 
0, f3g = 2 (Davis et al. 1985). Dashed curves: Qg = 8.4 x 10^ erg g-^ f3b = 0, (3g = 0.5 
(Davis et al. 1985). The light lines have Qb = 10^ erg g^^; the heavy lines have Qb = 10^ 
erg g~^. The horizontal lines indicate the collision energy for KBOs with e = 0.001 (light 
grey), KBOs with e = 0.01 (medium grey), KBOs with e = 0.1 (heavy grey). 
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Fig. 2. — Removal timescale as a function of size for KBO collisions with e = 0.04. The 
legend lists log Qy, and (ig for each curve. 
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Fig. 3. — Break radius, r?,, as a function of the bulk strength, log Qb, for collisions with e = 
0.04. The legend lists Pg for each curve. Heavy lines have Cg = 2.25 x lO^(^-'^^^i^s) erg g~^; 
light lines have Cg = 0.225 x lO^(^-'^^~Ps) erg g~^. The upper set of curves plots results for a 
minimum mass solar nebula; the lower set plots results for models with 1% of the mass in a 
minimum mass solar nebula. 
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Fig. 4. — As in Figure 4, for e = 0.2. 
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Fig. 5. — Size distributions at 4.5 Gyr for numerical models with constant e = 0.04. The 
legend lists log Qb for each model. 
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Fig. 6. — Variation of r^ with Qb for numerical calculations at 40-47 AU with constant 
e=0.04. The legend indicates the initial mass in solids for each set of calculations. 
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Fig. 7. — Size distributions at 4.5 Gyr for complete numerical models of KBO evolution at 
40-47 AU. The legend lists log Qb for each model. 
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Fig. 8. — Variation of r^ with Qb for complete numerical calculations of KBO evolution at 
40-47 AU. All calculation begin with a mass in solids equivalent ot the minimum mass solar 
nebula. The legend indicates log Qb for each set of models. Calculations with Neptune use 
a simple model for the growth of Neptune at 30 AU. 
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Fig. 9. — As in Figure 7, for models with Neptune at 30 AU. 



-34- 



w 

c 

0) 
Q 

E 

3 



2 - 



-2 









-4 - 



-6 - 



^ -8 



-10 - 



-12 - 




•''''^^ 







W. 



D 



o 



^Vv, 



Wvl 



D 



• 3.0, 0.0, 0.15, 1.25, N 

V 7.3,-0.4,1.5,1.25 

O 3.0,0.0,1.5,1.25 

D 3.0,0.0,0.15,1.25 






J. 



J. 



J. 



± 



20 



25 



30 35 

R Magnitude 



40 



45 



Fig. 10. — KBO luminosity functions derived from the planet formation model. The legend 
indicates fragmentation models for each model; the model with 'N' has Neptune stirring. 
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Fig. 11. — As in Figure 10, with observations of KBOs added for comparison (Bernstein et 
al. 2003). 
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Fig. 12. — Fraction of initial mass remaining in complete KBO models at 4.5 Gyr. Filled 
circles: models without Neptune. Open circles: models with Neptune at 30 AU. At the 
bottom of the plot, the horizontal dashed line indicates the ratio of the current mass in the 
Kuiper Belt to the mass in solids of a minimum mass solar nebula. 



